EFFICIENT AND LONG-TIME ACCURATE SECOND-ORDER 
METHODS FOR STOKES-DARCY SYSTEMS 
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Abstract. We propose and study two second-order in time implicit-explicit (IMEX) methods for 
the coupled Stokes-Darcy system that governs flows in karst aquifers. The first is a combination of a 
second-order backward differentiation formula and the second-order Gear's extrapolation approach. 
The second is a combination of the second-order Adams-Moulton and second-order Adams-Bashforth 
methods. Both algorithms only require the solution of two decoupled problems at each time step, one 
Stokes and the other Darcy. Hence, these schemes are very efficient and can be easily implemented 
using legacy codes. We establish the unconditional and uniform in time stability for both schemes. 
The uniform in time stability leads to uniform in time control of the error which is highly desirable for 
modeling physical processes, e.g., contaminant sequestration and release, that occur over very long 
time scales. Error estimates for fully-discrctized schemes using finite element spatial discretizations 
are derived. Numerical examples are provided that illustrate the accuracy, efficiency, and long-time 
stability of the two schemes. 
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1. Introduction. Karst is a common type of landscape formed by the disso- 
lution of layers of soluble bedrock, usually including carbonate rock, limestone, and 
dolomite. Karst regions often contain karst aquifers which are important sources of 
potable water. For example, about 90% of the fresh water used in Florida comes 
from karst aquifers |27| . Clearly, the study of karst aquifers is of great importance, 
especially because they are seriously threatened by contamination [29] . 

A karst aquifer, in addition to a porous limestone or dolomite matrix, typically 
has large cavernous conduits that are known to have great impact on groundwater flow 
and contaminant transport within the aquifer. During high-rain seasons, the water 
pressure in the conduits is larger than that in the ambient matrix so that conduit- 
borne contaminants can be driven into the matrix. During dry seasons, the pressure 
differential reverses and contaminants long sequestered in the matrix can be released 
into the free flow in the conduits and exit through, e.g., springs and wells, into surface 
water systems. Therefore, the understanding of the interaction between the free flow 
in the conduits and the Darcy flow in the matrix is crucial to the study of groundwater 
flows and contaminant transport in karst region. 

The mathematical study of flows in karst aquifers is a well-known challenge due 
to the coupling of the flow in the conduits and the flow in the surrounding matrix, 
the complex geometry of the network of conduits, the vastly disparate spatial and 
temporal scales, the strong heterogeneity of the physical parameters, and the huge 
associated uncertainties in the data. Even for a small, lab-size conceptual model 
with only one conduit (pipe) imbedded in a homogenous porous media (matrix), 
significant mathematically rigorous progress has only recently been achieved. For the 
coupled Stokes-Darcy model that includes the classical Beavers- Joseph [5] matrix- 
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conduit interface boundary condition, see [7J [HI [H] • For various simplified interface 
conditions, see, e.g., [711111130]. Nonlinear interface conditions have also been proposed 
for Navier-Stokes/Darcy modeling; see, e.g., [TT1 118j . 

Due to the practical importance of the problem of flow and contaminant transport 
in karst aquifers, there has been a lot of attention recently paid to the development 
of accurate and efficient numerical methods for the coupled Stokcs-Darcy system; 
see, e.g., [TU1 EH [301 (30] among many others. The efficiency of the algorithms is a 
particularly important issue due to the large scale of field applications. Because of 
the disparity of governing equations and physics in the conduit and matrix, domain 
decomposition methods (also called partitioned methods by some authors) that only 
requires separate Stokes and Darcy solves seems natural; see, e.g., fl~2" l fT3 1 H6 l H7 1 l28 l 
I3T1 1321 1351 [37] . On the other hand, long-time accuracy of the schemes is also highly 
desirable because the physical phenomena of retention and release of contaminants 
takes place over a very long time scale. Therefore, there is a need to ensure the long- 
time accuracy of the discretization algorithms in addition to the standard notion of 
accuracy on an order one time scale. 

The purpose of this work is to propose and investigate two types of numerical 
methods for the coupled Stokcs-Darcy system. We discretize the system in time via 
either a combination of second-order BDF and and Gear extrapolation methods or a 
combination of second-order Adams-Moulton and Adams-Bashforth methods. These 
algorithms are special cases of the implicit-explicit (IMEX) class of schemes. The 
coupling terms in the interface conditions are treated explicitly in our algorithm so 
that only two decoupled problems (one Stokes and one Darcy) are solved at each time 
step. Therefore, these schemes can be implemented very efficiently and, in particular, 
legacy codes for each of the two components can be utilized. Moreover, we show 
that our schemes are unconditionally stable and long-time stable in the sense that 
the solutions remain uniformly bounded in time. The uniform in time bound of the 
solution further leads to uniform in time error estimates. This is a highly desirable 
feature because one would want to have reliable numerical results over the long-time 
scale of contaminant sequestration and release. Uniform in time error estimates for 
fully discrete schemes using finite element spatial discretizations are also presented. 
Our numerical experiments illustrate our analytical results. 

Our work can be viewed as a time-dependent non-iterative version of the steady- 
state domain decomposition work in |13[ 116] and as a generalization of the first-order 
schemes in [101 137] that achieve the desirable second-order accuracy withtout increas- 
ing the complexity. The backward differentiation-based algorithm can be viewed as an 
infinite-dimensional version of the scheme presented in [33] , but with the additional 
important result on time-uniform error estimates. The Adams- Moulton/Bashford 
based algorithm is new so far as the Stokes-Darcy problem is concerned. To the best 
of our knowledge, our uniform in time error estimates are the first of their kind for 
Stokes-Darcy and related systems. 

The rest of the paper is organized as follows. In Sj2j we introduce the coupled 
Stokes-Darcy system and the associated weak formulation as well as the two second- 
order in time schemes. The unconditional and long-time stability with respect to 
the L 2 norm are presented in $3j Section [4] is devoted to the stability with respect 
to the H 1 norm. The H 1 estimates are important for the finite element analysis; 
this is another new feature of our work, even for first-order schemes. In §5\ we 
focus on the error analysis of the fully discretized scheme using finite element spatial 
discretizations. Numerical results that illustrate the accuracy, efficiency, and long- 
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time stability of our our algorithms are given in 
concluding remarks in [J7] 



We close by providing some 



2. The Stokes-Darcy system and two types of IMEX methods. 

2.1. The Stokes-Darcy system. For simplicity, we consider a conceptual do- 
main for a karst aquifer that consists of a porous media (matrix), denoted by fl p £ M. d , 
and a conduit, denoted by ilf £ M. d , where d — 2, 3 denotes the spatial dimension. The 
interface between the matrix and the conduit is denoted T. The remaining parts of 
the boundaries of VL p and f2/ are denoted by dVL p and dttf, respectively. See Fig. 12.11 



-3a- 



Fig. 2.1. The physical domain consisting of a porous media Q p and a free-flow conduit Qf. 

The coupled Stokes-Darcy system governing fluid flow in the karst system is given 
by [HUB] 



S-£ - V • (KV0) = / 



at 



-V-T(u/,p) = f and V-u /= in £l f , 



(2.1) 



where the unknowns are the fluid velocity Uf and the kinematic pressure p in the 
conduit and the hydraulic head <f> in the matrix; the velocity in the matrix is recovered 
from Up = — KV0. In (|2.1j) . f and / denote external body forces acting on the domains 
Clj and rip respectively, and T(v,p) = v(Vv + V T v) — pi denotes the stress tensor. 
The parameters appearing in (|2.1|) are the water storage coefficient S, the hydraulic 
conductivity tensor K, and the kinematic viscosity of the fluid v. 

For simplicity, we assume homogeneous Dirichlct boundary conditions for the 
hydraulic head and fluid velocity on the outer boundaries d£l p and d£l / , respectively. 
On the interface T, we impose the Beavers-Joseph-Saffman-Jones interface conditions 



- T J ■ ( T ( U /,P/) ' n /) = OLBJSjTj ■ Uf, j = 1, 

-ri/ • (T(u f ,p f )-n f )=g<j>, 



(2.2) 



where nj denotes the outer unit normal vector for flf and {T ? }j = i ; 2,...,d-i, denotes 
a linearly independent set of vectors tangent to the interface T. The additional pa- 
rameters appearing in (|2.2p are the gravitational constant g and the Beavers-Joscph- 
Saffman- Jones coefficient cxbjsj- 
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2.2. Weak formulation. We denote by an( l II ' \\d the standard L 2 (D) 
inner product and norm, respectively, where D may be Qf, Sl p , or T. We often 
suppress the subscript D if there is no possibility of confusion. We define the spaces 

H/ = {ve (H\n f )) d | v = oonan / \r} 

H p = {i> G H\n p ) I iP = ondn p \T} 
Q = L 2 (n f ), W = H/X H p . 

Dual spaces are denoted by (•)' and duality parings between spaces and their duals 
induced by the L 2 inner product on the appropriate domain are denoted by (•,•). 

A weak formulation of the Stokes-Darcy system (|2.1I) is derived by multiplying 
the three equations in that system by test functions v G Hj, gip G H p , and q G 
Q, respectively, then integrating over the corresponding domains, then integrating 
by parts the terms involving second-derivative operators, and then substituting the 
interface conditions (|2.2[) in the appropriate terms. The resulting weak formulation is 
given as follows [7J[T5]: given / G (H p )' and F G (H/)', seek G H p , u/ G H/, and 
p G Q, with d(j)/dt G (Hp)' and du/dt G (H/)', satisfying 

((u t , v)) + a(u, v) + 6(v,p) + o r (iI, v) = («f , v))) 

6(u,g) = 0, 

where u = [u, 0] T , v = [v,-0] T , and f = [i,gf] T and where (-) t = d(-)/dt. In (|2.3|) . 
we have that 

((ut,v)) = (u t ,v) 0/ +gS{(f>t,tp)n p , b(v,q) = -(q, V ■ v) J2/ 

a(u, v) = a/(u, v) + a p (0, V) + aBJS,/(u, v) (2.4) 

a r (u,v) =5(0,v-n / ) r -.9(u-n / ,V')r, (((f,v))) = (f,v) n/ + g{f,^)n p , 
where 

a / ( u > v ) = f(Vu, Vv) sl/ , Op(<j>, ip) = 5(KV0, V-0)o p 
asjs,/(u, v) = a B j(u ■ f , v • f) r . 

In (|2.3p . Uf, 0, and p are the primary variables; as mentioned before, once the hy- 
draulic head <f> is known, one can recover u p , the velocity in the porous media, via the 
Darcy relation u p = — KV</>. 

It can be shown that the bilinear form «(•,•) is coercive; indeed, we have that [7] 

a(u,u) > (j/||Vu|| 2 + .gX min ||V0|| 2 + a BJ ||u-f||2) >C a ||Vu|| 2 , (2.5) 

where C a = min(f, gK min) > an< i where K m - m denotes the smallest eigenvalue of IK. 
We define the norms ||u|| a = a(u, u)a and ||v||s = ((v, v)) 5 . We have that ||v||s is 
equivalent to the L 2 norm, i.e., we have that 

c.||v||s < ||v|| <Cs||v|| s . (2.6) 

2.3. The second-order backward-differentiation scheme (BDF2). The 

first scheme we introduce discrctizes in time via a second-order BDF whereas the 
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interface term is treated via a second-order explicit Gear's extrapolation formula. We 
propose the following algorithm: for any v £ W and q £ Q, 



(2.7) 



(( 3fl " +1 2At + — >*)) + a ^ l+1 ^ + b (^P n+1 ) + v) 

= (((f"+ 1 ,v)))-S r (2u n -u"- 1 ! v) 
b(u n+ \q) = 0, 

where the artificial stabilizing term a s t(-, ■) is defined as 

a st (u,v) =7 / (u-n / ,v-n/)r+7p(^,V , )r (2.8) 

with parameters 7/,7 P > and ar(u, v) is defined as 

a r (u, v) = a r (u, v) - a st (u, v) . 

2.4. The second-order Adams-Moulton-Bashforth method (AMB2). For 

the second scheme, we combine the second-order implicit Adams-Moulton treatment 
of the symmetric terms and the second-order explicit Adams-Bashforth treatment of 
the interface term to propose the following second-order scheme: for any v £ W and 
q £ Q, 



u" 



, v) ) + a (D a u n+1 , v) + b (v, D aP n+1 ) + a st (D a u n+1 , v) 



At 

= (((f™+i,S?)))-a r (|u 
b(D a u n+ \q) =0 



where D a denotes the difference operator that depends on a parameter a and is defined 
by D a v n+1 = av n+1 + (§ - 2a) u" + (a- |) w"" 1 . The stabilizing term a st (-, •) is 
defined as in 



2.5. Efficiency of the schemes. The implemented schemes are highly efficient 
because we can decouple the Stokes and Darcy subproblems: 

1. given u™, u n_1 

2. set v = [v, 0] so that all terms involving <j) drop out and we only need to use 
a fast Stokes solver to obtain u n+1 ; 

3. set v = [0, if)] so that all terms involving u drop out and we only need a fast 
Poisson solver to obtain cf> n+1 ; 

4. Set n = n + 1 and return to step 1. 

Note that steps 2 and 3 can be solved independently. Moreover, legacy codes can be 
used in each of those steps. 

3. Unconditional and long-time stability. The goal of this section is to 
demonstrate the unconditional and long-time stability, with respect to the L 2 norm, 
of the two second-order schemes proposed in ^2j We first recall a few basic facts and 
notations that are needed below. 

Recall that the G-matrix associated with the classical second-order BDF is given 

by 



G 
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with the associated G-norm given by ||wj| G = (w, Gw) Vw 6 (L 2 (f2)) 2 . The follow- 
ing identity is well-known (see, e.g., [13]): for any Vi £ L 2 (il), i = 0, 1, 2, 

3 „ , 1 \ 1 /,, ,|2 || ||2 1 | 11^2 - 2v 1 +v \\ 2 
-v 2 - 2ui + -wo,«2 J = 2 (II w iIIg - II w o|Ig) + ^ > I 3 - 1 ) 

where Wo = [woj w i] T and wi = [«i,V2] T . We also apply the G matrix to functions 
belonging to W: for any w € W 2 , define |w| G = (w, Gw). Then, for any v, <G W, 
i = 0,1,2, 

3 - O-^ 1 - - \\ 1 l\ |2 I ,2 N , |V2-2V1 + V0|| 

-v 2 - 2vi + -v ,v 2 ^ = - (|wi| G - |w | G ) + 2., 

where w = [v ,-^ 1 ] 7 " and W! = [vi,v 2 ] T . 

The G-norms are equivalent norms on (L 2 (Vt)) 2 in the sense that there exists 
Ci,C u > such that 

G/||w|| G < ||w|| 2 <C u \\w\\% and Ci\\w\\% < |w| G < C u \\w\\ G . 

We next recall the following basic inequalities: 
• trace inequality: if v € W, then 



||v|| r < CtrVl|v||||Vv||, ||v|| r < G tr ||Vv|| (3.2) 

• Poincare inequality: if v € W, then ||v|| < Cp\\ Vv|| 

11 2^2 2 

• Young inequality: a^b^c < =|- + + V o, 6, c, e > 0. 
Other variants of Young's inequality will also be used. 

The following estimate follows from the basic inequalities. 

Lemma 3.1. Let a 7 (-, •) and a s t(-, ■) be defined as in (|2.4j) and (|2.8[) . respectively. 
Then, there exists a constant C c t such that 

|a st (u,v)| + |a r (u,v)| < G ct ||u|| r ||v||r Vu,v e W. 

Proof. By the definition (|2.8j) of a s t (■,■), we have 

|a st (u,v)| < 7/Ku-n/, v n / ) r | + J p \{4>, ip)r\ 

< 7/llu-n/Hrllv -n/Hr + 7 P ll0j|r||V'l|r (3.3) 

< 7ma X (||u • n/|| r ||v • n/llr + \\cj)\\ r \\ ip\\ r ) , 

where the triangle and Cauchy-Schwarz inequalities are used and 7 max = m a- x {7/) 7p}- 
Similarly, by the definition (|2.4p of etr(-, •), we have 

|o r (fi, v)| < g (||0||r||v ■ n/|| r + ||u • n/||r||^llr) . (3.4) 
Note that u = fu, 6] T so that 



u 



| 2 =||u|| 2 + ||0|| 2 = ||u.n / || 2 +||u.r|| 2 + ||0|| 2 



Then, combining (|3.3p and (|3.4p . we obtain 

|Oat(u,v)| + |o r (u,v)| 

< (7max(||u ■ n/|| r ||v • Hf\\r + U\\r\\ip\\r) + s(IM|r||v • n/|| r + ||u • n/||r||V'l|r)) 

< max{7 max ,g}(||u • n f \\ r + ||</>||r)(||v ■ n/|j r + H^llr) 

< v/2max{7 max ,g}||u|| r ||v|| r . 
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The lemma is proved by setting C ct = y/2 max{7 max , g}. □ 

For the sake of brevity, we introduce the BDF difference operator Dv n+1 = 
l^n+i _ 2 v n _|_ iij"^ 1 and the central difference operator Sv n+1 = v n+1 — 2v n + 

3.1. Unconditional stability of the the BDF2 and AMB2 schemes. 

3.1.1. Unconditional stability of the BDF2 scheme. 

Theorem 3.2. Let T > be any fixed time. Then, the BDF2 scheme (|2.7|) is 
unconditionally stable on (0,T]. 

Proof. Setting v = u n+1 = (u n+1 , (f> n+1 ) in the BDF2 scheme pTTj) . we have 

i((flu"+\u" +1 ))+ B (if+ 1 ,u" +1 ) + Ssi (iu" +1 I u" +1 ) 

= (((f" +1 ,u" +1 ))) -a r (2u n -u"-\u" +1 ). 
From p. II) and the skew-symmetry of ar(-, we obtain 

l\^n\l-l\^n-l\ 2 G + 7II^U" + 1 ||| + Ato(u" +1 ,u" +1 ) + Ato it (u n+1 ,U n+1 ) 

(3.5) 



= At(^(((f™ +1 ,u n+1 ))) +5 r (-2u n + u"- 1 ,u" +1 ) 

where w„ = [u" +1 ,u"] T . Also, from the definition of the bilinear form a st (-,-), 
Lemma |3. 11 the trace inequality, and Young's inequality, we have 

a r {-2u n + u n -\u n+1 ) < C ct \\ -2u' 1 + u"- 1 || r ||u" +1 



r 

< C ct Cl\\ -2u" + u n - 1 ||*|| -2Vu n + Vu n - 1 ||*||Vu n+1 || 

- 2u" + u"- 1 i|^||Vu™ +1 || (V2||Vu n ||* + llVu"" 1 !^ 



< a*c 2 



<^|w„-i| 2 + ^!|Vu» +1 || 2 + ^||Vu'1 2 

+ ^|w«-il 2 + ^l|vd" +1 || 2 + ^livu"- 1 !! 2 . 

2 6 6 

For the forcing term, we have 

^3 M ?n-l-lii2 , C a 



(3.6) 



«<r+ 1 ,u"+ 1 )))<^||f"+ 1 || 2 + -^||u"+ 1 || 2 . (3.7) 



6C£ 



After we discard the nonnegative terms ||<5u' i+1 1|| and a s t(u' l+1 , u' l+1 ), noting that 
||Vu™ +1 || > ^||u" +1 ||, and using (J3HJ) and ([3770 . (j3~5|) becomes 



\w n \ 2 G + C Q Ai||Vu™ +1 || 2 < C 3 \\f n+1 \\ 2 At 

or 

+ (1 + (Ci + C 2 )At)|w„_ 1 | 2 ; + — ^=l||Vfi"|| 2 + :i ^||Vu 



2 . 2C a At 2 C a At Mr7 -j n _i M 2 



Next, by adding Cla 3 At || Vu"|| 2 to both sides of this inequality, we deduce 

C a At 2 (Ci+C 2 ) ., +1||2 C a M 2 {d +C 2 ) u vr n»2 

^ l+ (l + (C 1 + C 2 )At) i|Vu 11 + 3(l + (C 1 + C 2 )At) l|Vu 11 

< C 3 ||f n+1 || 2 Ai + (1 + (d + C 2 )At)E»-i 

< e (C 1+ C 2 )T So + C^_ e (C 1+ C 2 )T max|| fr l+ l ||2i 

Ci + G 2 " 
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where 



E n = |W„| G 



C a At 



Vu 



n+l II 2 



C a At 



IVu 



n||2 



(1 + (Ci + C 2 )At) 11 ' 11 3(1 + (d + C 2 )At) 1 

Thus, the unconditional stability of the BDF2 scheme is proved. □ 

3.1.2. Unconditional stability of the AMB2 scheme. We introduce the 
parameters 



(3.8) 









1 


ati = 


1- 




a-- 



fa = a x + a 2 , fa=2a- fa, fa = ^ (fa + fa)- 



Theorem 3.3. Let T > be any fixed time and let 1/2 < a < 1. Then, the 
AMB2 scheme (|2.9p is unconditionally stable in (0, T] . 

Proof. Setting v = u" +1 = (u' l+1 ,0™ +1 ) in pity , we deduce 



-, u ,l+1 )) + a(D a u n+ \ u™ +1 ) + a st (£> Q u n+1 , u' l+1 ) 



a st (-u --u ,u 



(3.9) 



Combining the two a S i(-, •) terms and using the basic equality 2(a — b)a = \a\ 2 — \b\ 2 

|2 



\a — b\ 2 , we have 
1 



At 



u 



™+ 1 l|2 _ || fjnil 2 , M.-tn + l _ -;n||2_ 



LS' 



- U 



u" +i - u"|||) + 2a (L> Q u" +i , u 



n+l fjn+n 



2(P+^,u" +1 ) - 2a r f-u" - -u"- 1 ,!! 



n+l 



2aa st (<5u" +1 ,u n+1 ). 



(3.10) 



Note that a > a±+a 2 provided | < a < 1. Therefore, /3i = 2a— /?3 = 2a— (ai+aa) > 
a and hence fa > fa > fa when ^ < a < 1. By the Cauchy-Schwarz inequality, we 
then have 

2a (D a u n+1 ,u n+1 ) > 2aa(u n+ \u n+1 ) - a 1 (a(u n+ \u n+1 ) + a{u n ,u n )) 

-a 2 (a(u n+1 ,u n+1 ) + a(u"- 1 ,u n - 1 )) (3.11) 
= faa(u n+ \ u' i+1 ) - a ia (u", u") - a 2 a(u"- 1 , u™" 1 ). 



Similarly as for p.6[) . for the interface coupling term, there exists a constant C5 such 
that 



2aa st (-2u" + u"- 1 ,u™+ 1 ) 



< 



Ca(^~ gg) || Vfln ||2 + 2^11^111 + g^LZ ^) ||Vfl"+l||2 

+ g^L^ llVu"- 1 !! 2 + Csllu-- 1 !! 2 , + ^_^||Vu" +1 || 2 . 



(3.12) 



For the forcing term, there exists a constant Cg such that 

2A<(((f n +^u' 1+1 ))) < C 6 Ai||P+*|| 2 + 9^1f ^ 2) At||Vu" +1 || 2 . (3.13) 
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Substituting ([3"3l~ lt -([3TT5 |) into (j3TTU|) yields 

+ ^LzM At || V u"+ 1 |f 

,ls + 2aAta st (u" 

< C 6 ||r + 3|| 2 At + (1 + 2C 5 At)||u"|| 2 s + CsAtHu"- 1 !! 2 + AtoiHu"^ 

+ AtaaHu"- 1 ]! 2 + Ca ^\~ ^ At\\Vir\\ 2 + ° a{Pl ~ h) At|| Vu"' 1 1| 2 . 
4 8 

Define the energy 

TP - ll.r l ll 2 _|_ ^At |i rl -i||2 , feAt 2 Ata 2 || 1 - T n-l||2 

3gq(A -fe) A ,, vfln||2 Cq(/?i -fe) At || Va n-i|ia 
+ 8(l + 3C 5 At) 11 " + 8(l + 3C 5 At) l|V 11 ' 

Then, discarding the last two positive terms on the left-hand side of (|3. 14[) and adding 
C 5 Ai||u n ||| + a 2 At|ju"|j 2 + ^(ft-fe) At||Vu"|| 2 to both sides, we obtain 

3C B At a a (/3 2 -/? 3 )At + 3C 5 feAt 2 +1 a 

En+1 + l + 3C 5 At l|u lls + l + 3C 5 Ai I|U L 

, 3C 5 a 2 At 2 a ZCsCajh - fa)At 2 2 



l + 3C* 5 Ai" " a 8(l + 3C 5 At) 
C a (/3i-/3 2 )At + 12C 5 C a (^ - /3 2 )At 
8(1 + 3C 5 At) 

< CellP+s || 2 At + (1 + 3C 5 At)^ 



2 



|Vu 



n+l||2 



Discarding all terms on the left-hand side, all of which are positive, except for E n+ \, 
we are left with 

E n+1 < C 6 Atmax||f n+ 5|| 2 + (1 + 3C 5 Ai).E n . 

n 

Then, by recursion, 

E n < e^ T E 1 + -g-e 3C5T max||r+^|| 2 

3C.5 i 

so that the proof of the theorem is complete. □ 

3.2. Long-time stability of the the BDF2 and AMB2 schemes. 

3.2.1. Uniform in time estimates for the BDF2 scheme. 

Theorem 3.4. Assume that f G L°°(L 2 (tt)) and that the time-step restriction 
(I3.19[) is satisfied. Then, the solution to the BDF2 scheme (|2.7p is uniformly bounded 
for all time. Specifically, there exist < Ai < 1, A 2 < oo, and Eq > such that 

\\u n \\ 2 < C U A?S + A 2 . 

Proof. Recall that o r (u n+1 , u n+1 ) = 0. Therefore, by LemmaEj] 
a r (-2u" + u"- 1 ,u" +1 ) - a st (<5u" +1 , u n+1 ) 

= a r (6u n+ \u n+1 ) < c ct ||<5u' l+1 || r ||u" +1 || r - ^ 
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The trace and Poincare inequalities imply 



|5u ,l+1 || r ||u" +1 || r < Cl\\6u n+1 \\^\\\75u n+1 \\i\\\7u n+1 \ 



(3.16) 



< CjC 2 J<5u n+1 |||(||Vu n+1 ||^ + V2||Vu"||' + ||Vu ri - 1 ||^)||Vu ,l+1 || 
The three terms on the right-hand side can be bounded using Young's inequalities: 



|<5u" +1 |||j|Vu" +1 ||s < |||Vu n+1 || 2 + p||M' l+1 i|| 



V2\\6u 



■n+l|| 2 \l\7ii n 



Vu"|| 5 l|Vu n+1 || < -||Vu 
8 



;n+l||2 



A!|Vu"|| 2 + ^yu" +1 || 2 



16 



|<5u"+ 1 |||||Vu"- 1 ||2||Vu™ +1 || < |||Vu" +1 || 2 + ^UVu"- 1 !! 2 + ^||<5u ,l+1 |||. 



Then, setting s 
(f3TT6|) that 



-o — - r 



Cs C c tC + . 



-, we deduce from these three inequalities, (|3.15(1 . and 



a r (5u n +\u n+1 ) < &|Vu n+1 | 



C a 



|Vu"| 



16 11 11 16 
The forcing term can be bounded as 



a 



(3.17) 



2Ci 



«<r+ 1 ,u«+ 1 )))<^||f»+ 1 || 2 + 



Cg 

C a " 8C?> 



?n+l 



(3.18) 



Combining (|3~5l> and ([H} with ff3TT7|> and (f3~T8]) . we obtain 



|w„| 2 +C a At||Vu n+1 | 



1 268C s 2 C c 4 t C t 8 . 

2 C3 



<5u 



n+l I 



< 



4C 2 Ai, 



f" +1 || 2 + |w„-i| 2 G + ^ ||Vu"|| 2 + ^ IIViT 



If the time-step restriction 



At < 



C;: 



(3.19) 



is satisfied, this leads to 



|w n | 2 3 + G I At||Vu" +1 || 2 



< 



4C 2 At, 



C a 



r +1 f + \*n-i\ 2 G + ^ iivu'i 2 + ^ llvu-i 2 



Adding 3C g At ||Vu"|| 2 to both sides of the above inequality, we obtain 



|w„| 2 ; + C a At||Vu" +1 || 2 + ^^ 



< 



4 C 2 At 

C a ' 



JVu" 
C a At 



r+T + l^-il 2 G + ^llvfii 2 



C a At 

8 



Ivu"- 1 !! 2 
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which is equivalent to 

E n + ^-At ||Vu" +1 || 2 + ^Ai ||Vu"|| 2 < E n . x + l^Hr+iH 2 , (3.20) 

where E n = \tf n \% + %^ ||Vu" +1 || 2 + ||Vu"|| 2 . 

Utilizing the Poincare inequality and the equivalence of the G-norm and the L 2 - 
norm, we have 

^ ||Vu"+i 2 + °± UVu-ll 2 > °± ||Vu"+i 2 + ^ ||Vul 2 + ^Iwnl 2 ,- 
Therefore, setting C-j = min{ C ^i a , 53*}' wc navc lrom <|3-20p that 
(1 + C 7 At)E n < B n _! + l^||f«+ 1 || 2 . 



A simple induction argument leads to 



1 \" , 4G 2 (l + G 7 At) f 2 

S„ < — — E n H „ „ max P r. 



l + G 7 Ai7 " C a G 7 
Recall that ||u"|| < C u E n . Hence, the theorem is proved with Ai = y+U~M anc ^ 

The following corollary is used in the analysis of the fully-discrete BDF2 scheme; 
see gSU 

Corollary 3.5. In addition to the assumptions of Theorem 13.41 assume that 
the second time-step restriction (|3.22[) is satisfied. Then, 

||u™|| 2 < CA^ 2 (||u°|| 2 + ||u 1 || 2 + At 2 ||Vu°|| 2 + At 2 ||Vu 1 || 2 ) +CA 2 . (3.21) 

Proof. For the interface term (|3.15|) . we can derive another estimate. From 
(|3~^|) and noting that ||V5u" +1 ||5 < V2(|| Vu n+1 j| 3 + ||Vu"||i + HVu"- 1 ^) and 
||u Tl+1 ||2 < v/2(||5u" +1 ||4 + ||u n ||5 + Hu"- 1 !!!), we have 



a r (<5u n+1 ,u n+1 ; 



n+l 



< CIIM'^HI J2 HVff ||^(|| ( 5u"+ 1 ||| + ||u"||J + ||u"- 1 ||I)||Vu"+ 1 ||i 

3=n— 1 
n+l 

= C\\Su n+1 \\ s l|Vu J ||*||Vu" +1 ||^ 

j=n—l 
n+l 

+ c\\su n+1 \\l J2 IIViP'||*(||ir||| + ||tr- 1 |||)||Vir +1 ||i -.= s 1 + s 2 , 



j-n-l 

^2 

inequalities: 

n+l 



where C = 2C s C c tC 2 r . The terms in the right-hand side can be bounded by Young's 



^ / I ~ \ 

Si< E U^ll^ n+1 || 2 5 + 3C 2 At(||Vu" +1 || 2 + ||Vtf|| 2 ) 

Su n+1 \\ 2 s + S^AidlVu 11 " 1 !! 2 + ||Vu"|| 2 + 2||Vu n+1 || 2 ) 



j=n-l 
J-| 

8 At 1 
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and 



71+ 1 71 



s 2 < E E (^ii^ +1 iil + §iiva" +1 ii 2 +^iiu fc iii + ^iivff'ii 2 ) 

j=n—l k=n—l 

-i n 97/^2 n +! 

-^ll^Hl + ^llVu'^f + E ^"11**111 + E C 2 At||Vtff. 

fc= n— 1 j—n—1 



Now, if we require 



the interface term can then be bounded by 

S r (M" +1 ,ii" +1 )< i ^||5u" +1 ||| 

97C' 2 ~ 
+ Yl|Vti" +1 || 2 + — (IKHl + llff-li) + 4C 2 At(||Vu"|| 2 + UVu- 1 !! 2 ) 

which leads to another recursion formula: 

iw»i 2 G + ^ iivu" +i ir < ^^iif" +i ii 2 

+ \*n-i\ 2 G + ^77^(l|u n || 2 s + llu"- 1 !! 2 ) + 8C7 2 Ai 2 (||Vu"!| 2 + UVu- 1 !! 2 ). 

<~'a 

Using this relationship, it is easy to verify that 

E n < CAt(\\i n+1 \\ 2 + ||P|| 2 ) + C\V n - 2 \% + CA^dlVu"" 1 !! 2 + !|Vu"- 2 || 2 ). 
Specifically, for n = 2, we have 

E 2 < CAt(||f 2 || 2 + ||f 1 || 2 ) + qK,| 2 , + CA^ 2 (||Vu || 2 + ||Vu 1 || 2 ). (3.23) 

Combing with Theorem 13.41 completes the proof. □ 

3.2.2. Uniform in time estimates for the AMB2 scheme. We start with 
the following estimate. 
Lemma 3.6. Let 

E T = -2a r Qu" - ^u n -\xi n +A - 2aa st (Su n+ \u n+1 ) . 
Then, with B\ and B 2 defined in (|3.8[) , we have the bound 

^ < 4C^-A) || V3 n + l||2 + 2C (A-fe) ||vfln||2 

C (8 -B ) 9 (3 - 24) 

. Wl.Pl Pi) n _ ||2 * tri \ r< Ml.-tn+l .-ini|2 , r,^ ll.-sn .-*n-l||2 



9 



IVu™- 1 !^ + (C 8 + C 9 )||u™ +1 - u"||£ + 2C 9 ||u n - u"" 1 !^. 



Long-time accurate 2nd-order methods for Stokcs-Darcy systems 



13 



Proof. Recall that cir(")") is skew-symmetric. Therefore, 

\£ r \ < \2a r (u ,l+1 - u", u rl+1 ) - a r (u™ - u"" 1 , u" +1 ) | 

+ | - 2aa st (u™ +1 - u™, u" +1 ) + 2aa st (u n - u" -1 , u' l+ 

< 2C ct \\u n+1 - u' l ||r||u" +1 ||r + 2C ct \\u n - u"- 1 ||r||fl"+ 1 || r 

< 2C ct Cl\\u n+1 - u"|| 1/2 ||V (u" +1 - u") || 1/2 ||Vu n+1 || 
+ 2C ct Cl\\u n - u™- 1 !! 1 / 2 !^ (u" - if" 1 ) || 1/2 ||Vu'" +1 j| 

< 2C ct Cl\\u n+1 - u"|| 1 / 2 ||Vu n+1 || (||Vu" +1 || 1/2 + ||Vu"|| 1/2 ) 

+ 2C ct C 2 r ||u"-u"- 1 || 1/2 ||Vu' l+1 || (||Vu"|| 1/2 + UVu"- 1 !! 1 / 2 
Then, by Young's inequality and (|2.6[) . 

\£r\ < Ca(/?1 ~ ^l ||Vfl" +1 || 2 + C s \\u n+1 - u"||| 
9 

+ Ca(/3l-/32) ||vfln+1||2 + g «(ft-ft) || Vjr ||2 + C^+l , 



(3.25) 



9 

C«(ft - ft) „ Vir+ i|| 2 + -»^ -^^ ||Vu"|| 2 + C 9 ||u" - u'- 1 !! 2 

9 

V3 n+l||2 + Ca(/3l-ft) ||vfln _ 1||2 + ^ 



9 

C (ft " ft) 



9 

4C a (ft-ft) +1 2 2C a (ft-ft) ||2 



C a (ft - 


ft) 


9 




C a (ft - 


ft) 


9 




a (ft - 


ft) 


9 




2C a (ft - 


ft) 


9 



9 

^ Vu"- 1 !! 2 + (C 8 + C 9 )\\u n+1 - u"|| 2 + 2C 9 ||u" - if 



C a (ft - ft) ... 



□ 

Theorem 3.7. Assume that 1/2 < a < 1, f e L°° (L 2 (Vt f)) , and that the time- 
step restriction (|3.28[) is satisfied. Then, the solution to the AMB2 scheme (|2.9[) is 
uniformly bounded for all time. Specifically, there exist < A3 < 1, A4 < 00, and 
.Ei > smc/i t/ia£ 

IK +1 || 2 <C S A?E 1 + A 4 . 



Proof. The interface term has been estimated in Lemma 13.61 The forcing term 
can be bounded as 

2(((r+2,u" +1 ))) < c 10 ||f"+^|| 2 + c °(^-^) ||vfl » +1||2 ^ (3 26) 

Combining (j3~TTj) . (|3~24)l . and (pOS]) . pTTO]) becomes 

||u" +1 || 2 + Aift||u' l+1 || 2 + ^AlM At||Vu" +1 || 2 
+ (l-(C 8 + C 9 )Ai)||u" +1 -u"|| 2 
< llu^lll + Ciolir+^lpAt + Atoillu'l^ + Atoallu' 1 - 1 !! 2 (3.27) 

2c a (ft - ft) Ai||Vt j, i||2 + a(ft-ft) Af||vfl „-i, |2 



2C 9 Ai||u n -if- 1 ^ 
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Now add CnAi||u n ||2 + C 12 c ° ( ^ 2) Af||Vu"|| 2 to both sides, require that 

/3 2 - ai > C n > a 2 , 2>Ci2>l, 
require the time-step restriction 

M < C^W (3 - 28) 



and set 

E n =||u»||| + At ( ai + C u ) \ml + (2 + C 12 )^fc^At||Vu"|| 2 

+ Atosllu"- 1 !! 2 + Ca{ ^~ ^^ AtHVu"- 1 !! 2 + 2C 9 At||u" - u"- 1 !] 2 . 
9 

Then, (pT27)) becomes 

E n+1 + At (ft - a x - C n ) + (2 - C 12 )^^l At|| Vu" +1 || 2 

+ At (C n - a 2 ) ||ul 2 + (C 12 - 1) Ca(/?1 Q ~ ^ At|| Vu"|| 2 ( 3 ' 29 ) 
+ (1 - (C 8 + 3C 9 )At)\\u n+1 - u"||| < C 10 ||f n+ 5|| 2 At + E n . 
Because there exists a constant C13 > such that 

C 13 ||fl" +1 ||5 < (2-C 12 )^i-^||Vu" +1 || 



' 18 




a\ - Cxi 


) 


- a 2 ) 




C-C81- 


02) 


18 




C a (Pi- 


ft) 



C 13 (2 + C 12 ) Caif3l Q fe) At 2 < (2 - d 

Ci3 a^_-ft) At2 ^ (Ci2 _ ~, Af> 

2C 13 C 9 At 2 < (1 - (C 8 + 3C 9 )At) 

we have from Q3.29[) that 

(1 + C 13 At)E n+1 < E n + dolir+^fAt. 

Thus, we have 

II^IIS < E n+1 < (-^-Ye! + ^od + C 13 At) ma ?t 



.l + C 13 At7 C w 

Ci (l+t 

l+C 13 At " 4 Ci, 



Setting A 3 = 1 and A 4 = C ' lo(1 + Cl3At) max^ \\f i+ ?\\ 2 , by (JUD the proof is 



complete. □ 

Remark 1. Similarly to Corollary 13. 51 m £/ie error analysis, E\ can be taken as 
CdwolG + A^dlVuolP + HVuiH 1 )) in TheoremMi 

4. H 1 ^!) stability of the schemes. The purpose of this section is to prove 
uniform in time bounds for the solutions to the schemes (|2.7j) and ()2.9p with respect 
to the if 1 (SI) norm. This additional estimate is needed for the estimation of finite 
clement element errors for for the fluid velocity and hydraulic head with respect to 
the iJ 1 (SI) and for the pressure with respect to the L 2 (S1/) norm; see £15.11 



Long-time accurate 2nd-order methods for Stokcs-Darcy systems 



15 



4.1. Uniform in time i7 1 (SI) bound of the BDF2 scheme. In this subsec- 
tion, we assume that the time-step restriction (|3.19[) holds. We introduce the notation 
Bu n+1 = ^(u n+1 -u"). 

Lemma 4.1. The first-order discrete time derivative of the BDF2 scheme (|2.7|) 
is uniformly bounded in time. Specifically, we have 

||du" +1 || 2 < CX'l + Cmax||dP|| 2 , (4.1) 



where the positive parameter X\ < 1 is defined in Theorem \3A[ 

Proof. For the BDF scheme (|2.7[) . we take the difference of the n and n + 1 level 
equations to obtain 

-L((Ddu n+1 , v» + a(<9u rl+1 ,v) + b(v,8p n+1 ) + a st (Sdu n+1 ,v) 
= ( ( (Bi n+1 , v)» + ar (~2du n + <9u"- 1 , v) . 
Now setting v = du n+1 and using the skew-symmetry of ar, we have 

-L((Ddu n+1 , Bu n+1 )) + a(Bu n+1 , Bu n+1 ) + a st (5du n+1 , Bu n+1 ) 
= (((Br +1 ,Bu n+1 ))) +a r (SBu n+1 ,Bu n+1 ). 

The rest proof is a verbatim copy of the proof of Theorem 13.41 with f replaced by Bi . 
□ 

A direct consequence of the Lemma 14.11 is the following result, once we realize 



that -£- t Du n+1 = §<9u n+1 - ±<9u™ and -^- t Su n+1 = Bu n+1 - 8u n . 

Corollary 4.2. Let u" be the solution to the BDF2 scheme (|2.7p . Then, 

||-Ldu" +1 || 2 + ll-L^u'^ 1 !! 2 < CA'^ + Cmaxliaril 2 . (4.2) 
At At i 

The following technical lemma is useful in deriving the uniform in time H l (£l) 
bound. 

Lemma 4.3. Let {a n } be a nonnegative sequence that satisfies 

a n +i < ciAt(a n + a„-i) + c 2 A" + c 3 for n = 1, 2, . . ., 
where Ci, i = 1,2,3, are positive numbers andO < A < 1. Moreover, if At < - — — , 

">>'•*' ' J (1 + V5)ci' 

then, 

a n+ i < T ^ r — + A" ( + a 1 + ^-±a») . (4.3) 

Proof. Define b n +i = a n+ i + ^ c\Ata n . Then, 

K+i < 1 + ^ c 1 Atb„ + c 2 A" + c 3 . 
A simple induction leads to 

b n+1 < £ {^^c.At)^ 1 { C2 y + c 3 ) + (i±^ Cl At)\. 
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Now if At < . 2 * , we have 

(1+V5)ci 

b n+1 < % + ^ + x n h. 

The desired bound on a„ +1 follows from this inequality, the definition of b n+ \ and &i, 
and the fact that CiAt < 1 under the assumption. □ 
Remark 2. If At < ^ , then (|4~3)) implies 

a n +i < 2c 3 + A n (2c 2 + ai + —— — a Q ). 



Theorem 4.4. The BDF2 scheme (|2.7[) is asymptotically stable with respect 
to the iJ 1 (SI) norm in the sense that the H 1 ^) norm of the solution is uniformly 
bounded in time. 

Proof. Set v = u" +1 in the BDF scheme (|2.7[) and use the skew-symmetry 
property of ar to obtain 

a(u n+1 , u n+1 ) = -^((L>u" +1 , u n+1 )) + a r {5u n+1 , u n+1 ) + (((f' l+1 , u" +1 ))). 

Note that 

-^ t ((Du n +\u n+1 )) + (((f"+\u" +1 ))) < C (j|^u» +1 || + ||f n+1 ||) ||Vu" +1 || 
and 

a r (Su n+ \u n+1 ) < ^C||^<5u' i+1 || + ^||V5u" +1 ||^ ||Vu' l+1 ||,. 

Using the coercivity condition (|2.5p and (|4.2p . we deduce 

||Vu ,l+1 || < C* 9 (\\-^Du n+1 \\ + \\-Ldu n+1 \\ + ||f ,l+1 || + At(||Vu™|| + ||Vu" 

< C 15 \f + Cis max(||9f i+1 || + ||r +1 ||) + C u At(\\ Vu"|| + UVu"- 1 !)), 

i 

provided that At < C a . The desired uniform in time estimate then follows from 

Lemma H3] with a n = ||Vu"||, cj = Cu, c 2 = C 1B , and A = A? = (1 + C 7 Ai)-3. 
Specifially, provided that the time step is small enough in the sense that At < 

\j (l + v / ^)C'l4 fViri f iTnn_cfnr\ PAnrlifiAn lvi T omrrnJ H ^lio Tmyifiorl i a ( l~t~ V^jC 



, the time-step condition in Lemma I4T3I is verified, i.e., 2X t 14 At < 



i. Hence by Lemma [ 

||Vu" +1 || < CXf +2C 15 max(||ar +1 || + ||r +1 ||). 



□ 
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4.2. Uniform in time H 1 (CI) bound for the AMB2 scheme. In this sub- 
section, wc assume that the time-step restriction (|3.28[) holds. 

Utilizing the same arguments as for Lemma 14. 11 wc can deduce that the discrete 
time derivative of the solution of (|2.9|) is uniformly bounded in time. 

Lemma 4.5. For the ABM2 scheme, we have 

||9u n+1 || 2 < CAJp 1 +Cmax||(9r+5|| 2 , 

i 

where the positive parameter A3 < 1 is defined in Theorem 13.71 
Lemma 4.6. Let a n be a nonnegative sequence and let 

a n +i < c 4 a n + c 5 a„_i + c^A' 1 " 1 + c 7 for n = 1, 2, . . ., 



where a, i = 4, . . . , 7, are positive real numbers and < A < 1. Let £1 = Yf£tM±L C4 



and £2 = V c 4+4 Cj _££ _ // C4 + C5 < 1 , then 



^ en, , e \ 1 c 6 (max(A,Ci)) n 1 c 7 

«n+i < £1 (01 + £ 2 a ) H 7- cX + r . (4.4) 

l-min(A,^ 1-a 

Proof. Note that £1 < 1 because C4 + C5 < 1. Letting 6„ + i = a n+ i + £20™ we have 
b n +i < tib n + c 6 A n_1 + c 7 . 

By induction, 



bn+l < + C 6 ^ + C7 & 

i=l i=0 

Because £1 < 1, then J27=o £1 — T^T an< ^ 

y^„-i A i < (max(A,Ci))"~ 1 
fc 1 "l-min(A^)' 

Now (|4.4p can be obtained. □ 

Theorem 3. The ABM2 scheme is asymptotically stable with respect to the 
H (CI) norm. 

Proof. From (j3~9)) . 



a(D a u n+ \u n+1 ) = (((r+Ku n+1 )))+£r - ((Bu n+ \u n+1 )), 
where £r is defined in Lemma 13.61 Note that 

a(A l u"+\u" +1 ) > /? 2 |K +1 !| 2 + a a C8i - /3 2 )||Vu" +1 || 2 - ax||u"|| 2 - a^X, 
and 

(((f' i +^,u" +1 )))-((au" +1 ,u"+ 1 )) < c (||5u ii+1 || 2 + \\f n+1 *\\ 2 )+ Ca{p \~ ^ ||va" +1 || 2 . 
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By following the proof in (|3.25[) with small modifications, we obtain 

n+l 

£ v < C , At*||Vu n+1 || \\du j \\ 1/2 (\\Vu j \\ 1/2 + UVu^ 1 !! 1 / 2 ) 

j=n 

< CJ2 \\&&'\\ a + Ca(/?1 ~ /?2)A ^ (4||Vu" +1 || 2 + 2||Vu"|| 2 + UVu"- 1 !! 2 ) 

j=n 



Define a ± = a x + 2{(ix ~ h) At 2 and a 2 = a 2 + {f3l ~ M At 2 . Then, 

/3 2 ||u» +1 || 2 < SxHu-ll 2 + Sallfi- 1 !! 2 + CAr 2 + Cmax(||ar^j| 2 + ||P + *|| 2 ). 
Now, if At is satisfies 



At < 



'3 082 -03) 



V ft -ft 

then /3 2 > 5i + 5^2 and then, by Lemma 14.61 



? n+ii|2 ^ t n/i|rfi||2 , tf _ii^0ii2^ , C*(max(A 3 , £i)) n 



ll^ii^ffdia^ii+^iiifiil 

/3 2 A 3 (l-min(£,£ 

Cmax t (||ar+^|| 2 + ||r + ^| 2 ) 
+ ^2(1-6) ' 

where 6 - ^ 5j+ gf 2+51 and \M±gpLgl . D 

5. Error analysis. In this section, we study the convergence of the fully-discrete 
BDF2 scheme, where spatial discretization is effected using finite element methods. 
A similar study yielding similar results can be done for the AMB2 scheme; however, 
for the sake of brevity, we omit that study. 

Let Hf t h C H/, H Pt h C Hp, and Qh C Q denote conforming finite element spaces. 
Let Wh = Hf.h x Hp,h- We assume that the mesh is regular and that the parameter 
h is a measure of the grid size. We use continuous piecewise polynomials of degrees 
k, k, and k — 1 for the spaces Hf,h, H Pi h, and Qh, respectively. See [14] for details 
concerning such finite element discretizations. We also assume that the fluid velocity 
and pressure spaces Hf.h and Qh satisfy the discrete inf-sup condition necessary for 
ensuring the stability of the finite element discretization; see [21] . 

Definition 5.1. The Stokes-Darcy projection Ph : W xQ — s- x Q h is defined 
as follows. For any u G W and p £ Q, let P/jiJ £ and PhP £ Qh denote the finite 
element solution of 

a(P h u,-v h ) + b(v h ,PhP) + a r {PhU,v h ) = a(u,v h ) + b(v h ,p) + a r (u,\r h ) 

b(P h u,q h ) = b(u,q h ) 

for all v/, £ Wh and q h £ Qh- 

It is easy to see that for any u £ W and p £ Q, the exist unique P/jU £ W/j 
and PhP £ Qh- Moreover, if we assume that u £ (H k+1 (Q,f)) d x H k+1 (Q p ) and 
p £ H k (n f ), then (see, e.g., [8]), 

||u - P h u\\ + h\\V(u - P h u)\\ < h k+1 (\\u\\ Hk+1 + \\p\\ H >). (5.1) 
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Remark 4. The estimate (|5.1[) and the optimal error estimates derived below 
assume that u G (H k+1 ) d x H k+1 which requires that the interface T be sufficiently 
smooth. In this case, the finite elements may need to be modified near the interface, 
e.g., by using isoparametric finite element approximations [14) . In any case, in this 
paper we assume that the optimal error order of convergence can be obtained for the 
steady-state Stokes-Darcy problems using the same grids and finite element spaces. 

5.1. Error analysis of the BDF2-FEM scheme. The fully discrete BDF2- 
FEM scheme is defined as follows: for n = 0, 1, 2 . . ., seek u^ +1 € and p^ +1 e Q h 
such that 

i((Du^\v fc >> + a(ur\^) + Kv h! pr 1 )+a st (u^ +1 ! v fe ) 

= <(<f" +1 ,v fc >>>-a r (2u^-i!r 1 J S? /l ) + a s4 (2u^-ur 1 ,v /l ) ( 5 - 2 ) 
b(u»+\q h )=0 

are satisfied for all v/,. <G W/j and <?/, e Qh- Note that for all v/, e W/, and g/, g Qh, 
the exact solution satisfies 

((u t ,v fc » +a(u,v fe ) + a r (u,v/ l ) + 6(v h ,p) = «(f,v h >» g ^ 

6(u,? fc ) = 0. 

Theorem 5.2. Assume that the solution of the Darcy-Stokes problem (|2.ip is 
sufficiently regular in the sense that 

u e ^(O^-H 1 ) n H 2 (Q,T;H k+1 ), 

that the time-step restrictions (|3.19p and (|3.22p are satisfied, and that the finite el- 
ement spaces are chosen so that the projection error bound (|5.ip holds. Then, the 
solution of the fully- discrete BDF2 scheme (|5.2p satisfies the error estimate 

\\u(t) u£|| 2 < \\P h u(t a ) u°|| 2 + \\P h u( tl ) ui\\ 2 + CAt(\\V(P h u(to) u?J!| 2 
+ \\V(P h u(h) - ui)\\ 2 ) + C(At 4 + h 2 ^). 

Moreover, if the solution of the Stokes-Darcy problem (|2.ip is long-time regular in the 
sense that 

u G ^"(O.oo^ 1 ) n W 2 '°°(0,oo;H k+1 ), 

then, there exists a constant C a and a generic constant C independent of At, h, or n 
such that the solution of the BDF2 scheme (|5.2p satisfies the uniform in time error 
estimates 

\\u(t n ) u£|| 2 < CXr 2 (\\P h n(t ) ul\\ 2 + \\P h u(h) ul\\ 2 ) 

+ CAt 2 Ar 2 (||V(P, 1 u(io) - u?J|| 2 + \\V(P h u(h) uDll 2 ) (5.4) 
+ C{At i + h 2{k+1 ^) Vn 



||v(u(t n+1 ) - ur 1 )!! 2 + \\ P {t n+1 )-pi +1 \\ 2 

< C\r 2 (\\dP h u(h) - Bui}] 2 + \\dP h u(t 2 ) - Bu 2 || 2 ) 
+ CAt 2 Ar 2 (||V(aP, l u(i 1 ) - MDf + \\V(dP h u(t 2 ) - 9u 2 )|| 2 ) 
+ C(At 2 + h 2k ) Vn 
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provided that At < C a . 

Proof. Let e™ = u(t„) — uJJ denote the error at the time t = t n . Then, from (|5.2I) 
and (|5.3[) . we have 

i-((^e-+ 1 ,v /l )) + ffl (e n+1 ,v h ) + b(v h ,p(t n+1 ) - pl +1 ) + a T (e n+1 ,v h ) 



(5e" +1 ,v^) = (« +1 ,v fe )) -o r (*a(tn+i).Vh) 



6(e"+\ 0/l ) =0, 



where = -u t (i„+i) + ^Du(* n+1 ). Let p" = u(i„) - P h u(t n ) and n = 

Phu(t n ) — uJJ. Then, 0™ E W/j and is discretely divergence free, i.e., 

6(0", g h ) = Vq h eQ h . (5.6) 

Because e" = 9 n + p" 1 , the error equation (|5.5j) can be recast as 

^ ( ( D8 n + 1 , v h ) ) + a(6 n + 1 , v h ) + a r (0 n+1 , v h ) - a r (5d n+ 1 , v„ ) 
= <« +1 ,v fe » - a r (<Su(t„ + i), v h ) - b(v h ,p(t n+1 ) - p% +1 ) 

- -^ t ((Dp n +\v h )) - a(fT +1 ,v h ) - a r (p" +1 ,v fc ) +Z r (Sp n+1 ,v h ) 
= ((^ +1 ,v h ))-a r (Su(t n+1 ),^ h )-b(v h ,P h p(t n+1 )-pl +1 ) 

-^(( J Dp ra+1 ,v h » + a r (5p re+1 ,v fc ). 

Setting v/,, = n+1 , noting that ar(0 n+1 ,6 n+1 ) = 0, and using (|5.6[) results in 



— ((D0 n+ \0 n+1 )) +a{0 n+ \0 n+1 ) -a r (5d n+1 ,0 n+1 ) = (« +1 ,0 n+1 )) 
-5 r (ou(t„ +1 ),0" +1 ) - ±-((Dp n +\O n+l )) +a r (Sp' l+ \e n+1 ). 



(5.7) 



Letting w„ = [0™+ 1 ,0™] T and £„ = |w„|| + ^||V0" +1 |j 2 + ^||V0"|j 2 and 
following the lines of the proof of Theorem 13.41 we have 



E n + ^fAt 



V6> 



n+1 



^At 
4 



V0' 



+ CAt(|K +1 || 2 + \\V6u(t n+1 )\\ 2 + \\—DfP+ l \\ 2 + ||Vop 
By Taylor's theorem with the integral form of the remainder, we have 



(5.8) 



^ +1 \\ 2 < CAt 3 / ||u ttt || 2 dt < CAt 4 ||u ttt || 2 



L°°(0,t„+i) 



(5.9) 



and 



||V<5u(t„ +1 )|| 2 < CAt 3 / j|Vu ttt || 2 dt < CAt 4 \\X7u ttt \\ 2 L . 



'(o,t» +I )- 



(5.10) 



Moreover, using ()5.1[) . we have 

H^^+ 1 H 2 <c^ +1 )||^u(t„ +1 ; 



< c- 



tfik+l) ftn+l 



At 



l|u t ||| fc+1 *<c/i 2 ( fe+1) ||u t || 



(5.11) 



tn-l 



L° c (0,t, 1+ i;ff fc + 1 ) 
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and 



\V6p 



Sn+l||2 



< Ch 2k \\Su(t n+1 )\\' z Hk+1 < Cti 2k At 



2k a +3 



< c/i 2fc At 4 ||a tt |||. 



Combining ([5T5 |) - (|5TT2"1) . we have 

n+l 



£ n + ^At£|v0* 



h 2(k+l) 



< Eo + ClAt 



{\\utt 



a (0,t n+ i;J?*+ 1 )- 

h||Vu ttt || 2 )dt 



(5.12) 



(5.13) 



\u t \\ z Hk+1 dt + h 2k At 4 



\u tt \\ 2 Hk+1 dt 



The desired finite time error estimate follows from this and the assumed bound on 
the projection error p™. 

For the uniform in time L 2 (VL) bound, we again use (|5.8|l - (|5.12|l to obtain 



E n + ^At 



V6» 



n+l 



2 n 
+ ^At 
4 



V0' 



< £„-i + CAi(At 4 ||u« t ||^ (0iOo) + At 4 1| va« t ni» (0iOO) 



(5.14) 



+ ^ 2(fc+1) ||Ut|lL«=(0,oo;ff fc +l) + ft - 2fc ^ <4 ||Utt|lL«=(0,oo;_ff'" + l) 

< S„_i + GAi(A^ 4 + /i 2(fc+1) ). 
Using the Poincare inequality and the definition of the G-norm, we have 



^At 
2 



2 C a 



> ^At 
~ 4 



V0 n+1 

" V 0n+1 



^At 

4 



V6T 



G„ 



Af 



V6T 



G 2 G a . 
"8G^' 



Then, with G7 defined as in Theorem 13.41 we have from (|5.14[) , 

(1 + C 7 At)E n < £„_i + CAt(At i + h^ k+1 ^) 
A simple induction argument then leads to 

1 \ n-2 



E n < 



. 1 + G 7 Ai 
< GA™' 2 (||0 1 || 2 



E 2 + C(At 4 + h^ k+1 ^ 



oD 11 2\ 



\o 



CA^AtdjV^ 1 !! 2 + ||V0°|| 2 ) + G(At 4 + h^ k+1 ^), 



where Ai is defined as in Theorem 13.41 The bound (|5.4p then follows from Corol- 
lary OEU 

The uniform in time iJ 1 (fi)-norm error estimate on the velocity and the L 2 (f2) 
error estimate on the pressure can be derived as well after we combine the technique 
used above with techniques from Indeed, from (|5.9p and (|5.10[) and using the 
triangle inequality, we have 



|0wj 



n+l||2 



< CAt 



f„+i 



\u ttt \\ 2 dt<CAt 2 \\u ttt \\ 2 Lo 



(0,t„ H 



(5.15) 
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and 



\VB5u(t n+1 )\\ 2 < CAt f " +1 ||Vu ttt || 2 ^ < GA^HViltttllioo^^). (5.16) 



(5.19) 



Moreover, by the definitions of Ph, d, and D, 

W-LdD^W* < Ch 2 ^\\±dDu(t n+1 )\\ 2 Hk+1 

h 2(k+l) /•«„+! (5.17) 

< jf W^AWdt < ch^\\u u \\l^ tn+i ^ k+1) 

and by the triangle inequality and (|5.12l) . 

||V^p" +1 || 2 < Ch 2k At^u tt \\l^ {0 tn+i . Hk+ly (5.18) 

We combine (|5.15|) - (|5.18[) with the stability proof of Lemma [4.11 with small modifi- 
cation for the initial steps; see Corollary 13.51 As a result, we obtain 

\\d6 n+1 \\ 2 < CA^GI^H 2 + WB&W 2 ) 

+ CAt 2 A^ 2 (||V<9F, l 1 || 2 + WVdPh&W 2 ) + C(At 2 + /i 2(fc+1) ). 
Note that 

||^M|| 2 < CAt a ||a tt ||i. (Pit()> ||^|| 2 < CAth^^l^^y (5.20) 

Combining (f5TT9|) and (fOO]) with (|5~9]) and ([5TTj) and following the proof of Theo- 
rem 031 we have from (|5.7[) 

|| V 0„+i|| 2 < C 7A I 1 l - 2 (||90 1 || 2 + ||d0*|| 2 ) 

+ CAt 2 \ , l~ 2 {\\\7dP h 8 1 \\ 2 + WVdPhfPf) + C(At 2 + h 2{k+1) ). 

After adding the estimate of || Vp" +1 1| (see (|5.1[0 . we obtain the bound for || V(u(f„+i) — 
n n h +1 )\\ 2 - 

The error estimate for the pressure — Ph\\ can be obtained by standard mixed 
finite element analyses; see [21]. □ 

Remark 5. The uniform in time estimates given above imply that the method 
can be used to obtain an approximate solution of the steady- state equations in case the 
forcing term is time independent. This follows because the truncation errors listed in 
()5.8|) - (|5.12j) vanish for the time-independent problem. Consequently, we have 

\\u',l u»|| 2 < CA?(|K - u°|| 2 + At\\V(n° h u°)|| 2 ) 

for the steady-state problem. 

In the steady-state case, the methods we study can be viewed as a domain de- 
composition method with the discrete time n playing the role of an iteration number; 
see |13j for a related scheme. The current scheme also enjoys an exponential rate of 
convergence as does the iterative scheme proposed in |13j . 

Remark 6. Note that the uniform in time error estimate for the velocity with 
respect to the H 1 ^) norm and the the uniform error estimate for the pressure are 
not second order in time. We do not know if this is an artifact of our approach. 
However, our numerical experiments in the next section suggest that the long-time 
convergence rate for the pressure approximation may very well be first order as the 
analysis suggests. 
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6. Numerical results. Using three numerical examples, we now illustrate the 
theoretical results of the previous section. 

As was done in previous work [T3] [37], we set y = (0,1) x (1,2) and il p = 
(0, 1) x (0, 1), with Clf and £l p separated by the interface T = (0, 1) x {1}. We choose 
the standard continuous piecewise-quadratic finite element space, defined with respect 
to the matrix domain fi p , for approximating the hydraulic head <$>. We also choose 
the Hood- Taylor element pair, defined with respect to the conduit domain fi/, i.e., 
continuous piecewise-quadratics and continuous piecewise-linear finite element spaces 
for the fluid velocity and pressure approximations, respectively. Uniform triangular 
meshes are created by first dividing the rectangular domains f2 p and Q f into identical 
small squares and then dividing each square into two triangles. For illustrating the 
short-time properties of our schemes, we set the final time to T = 1; for illustrating 
the long-time behavior, we set T = 100. 

We use three examples with exact solutions. Example 1 is taken from [37], Ex- 
ample 2 from [13] . and Example 3 is a slight modification of an example in [8]. To 
illustrate the accuracy of our schemes, we assume that the error is of the order 
0{h 9% + At 02 ). We set At = h e and quantify the numerically estimated order of 
convergence rg = min(#i, 662) with respect to h by calculating 



0(x, t) = [2 — tt sin(7ra;)] [1 — y — cos(7rj/)] cos t. 

The right-side data in the partial differential equations, initial conditions, and bound- 
ary conditions arc then chosen correspondingly. As done in |37j , we set the parameters 
7p = 7/ = g = S = v — otBj = 1 and K = I; also, we set a = 0.8 for the AMB2 
scheme. 

For Table 16.11 we set At = h and present results for several values of h; the 
results illustrate the second-order in time accuracy for </>, u/, and pf. We also notice 
that BDF2 has a significantly smaller error than AMB2, illustrating the advantage 
of the former over the latter scheme, at least for this example. For Tables 16.21 and 
16.31 At is chosen to be a power of h to illustrate the spatial convergence rates. The 
results in those tables indicate that the spatial accuracy seems higher than the third- 
order suggested by our analysis. The extra half order of accuracy may be attributed 
to super-convergence or super-approximation behaviors; see [12j for a study of this 
phenomenon for the steady state case. 

Example 2. We next set the exact solution to the steady-state solution [13] 



r e w log 2 



\\u2h,e " u exact \\i2 

\\uh,9 — U exact \\i2 




u f (x,t) 
0(x,t) 



— sin(27ry) cosx , 2H r-sin 2 (7ry) sinx 



7T L TT 







(e y — e y ) sinx. 
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e u 




h 


BDF2 


AMB2 


BDF2 


AMB2 


BDF2 


AMB2 


1/16 


5.76e-005 


3.43e-003 


8.26e-005 


l.llc-004 


1.15e-002 


4.11C-002 


1/32 


9.53c-006 


8.76e-004 


1.98e-005 


2.74e-005 


3.02e-003 


1.07e-002 


1/64 


2.35e-006 


2.21c-004 


4.85c-006 


6.79e-006 


7.73e-004 


2.71c-003 


1/128 


6.00c-007 


5.55e-005 


1.20c-006 


1.69e-006 


1.96e-004 


6.85e-004 


favg 


2.20 


1.98 


2.04 


2.01 


1.97 


1.97 



Table 6.1 

Relative error and order of accuracy with respect to the spatial grid size h for Example 1 at 
t = 1 and with At = h. 





e<f, 


e u 


Cp 


h 


BDF2 


AMB2 


BDF2 


AMB2 


BDF2 


AMB2 


1/8 


6.16e-004 


6.83e-004 


8.14e-005 


8.37c-005 


2.81c-002 


3.04c-002 


1/16 


5.39c-005 


6.46e-005 


7.67c-006 


7.86e-006 


7.71c-003 


7.93c-003 


1/32 


4.70e-006 


6.01c-006 


6.99c-007 


7.16e-007 


2.03e-003 


2.05e-003 


1/64 


4.13e-007 


5.51c-007 


6.26e-008 


6.41e-008 


5.22e-004 


5.24e-004 


r avg 


3.51 


3.43 


3.45 


3.45 


1.92 


1.95 



Table 6.2 

Same information as in Table\6l]but for At = ft 3 - 5 / 2 . 







e u 


Cp 


h 


BDF2 


AMB2 


BDF2 


AMB2 


BDF2 


AMB2 


1/8 


6.17e-004 


5.82e-004 


8.11c-005 


8.17c-005 


2.78e-002 


2.85c-002 


1/16 


5.40c-005 


5.21e-005 


7.66e-006 


7.69e-006 


7.65c-003 


7.73c-003 


1/32 


4.71e-006 


4.62e-006 


6.99c-007 


7.01c-007 


2.04e-003 


2.03e-003 


1/64 


4.13e-007 


4.09e-007 


6.26e-008 


6.28e-008 


5.22e-004 


5.22e-004 


favg 


3.52 


3.49 


3.45 


3.45 


1.91 


1.92 



Table 6.3 

Same information as in Table RTTl but for At = h? 



It is easy to see, after scrutinizing the error analysis, that third order instead of 
second order in time convergence is expected for <j) and Uf in this steady-state case. 
Our numerical results in Table 16.41 for which we have set At = h, suggest 3.5-order 
convergence for these variables. The convergence rates are consistent with the results 
of [13] for the steady-state case. Thus, it seems that the super-convergence results of 
[12] hold for our new temporal discretization schemes. The rigorous demonstration of 
the super-convergence rate can be accomplished by following the analyses of [12] and 
will be discussed in future work. 

Example 3. To illustrate the long-time behavior of our schemes, we use the following 
exact solution that is a slight modification of an example in [8] : 

u/(x,i) = {[x 2 y 2 + e~ y ], [~xy 3 + [2 - Trsin(Tra)]]) [2 + cos(27rt)] 

p/(x, t) = — [2 — 7rsin(7ra)] cos(27ry)[2 + cos(27rt)] 
0(x,i) = [2 - irsm(irx)][-y + cos(7r(l - y))][2 + cos(27rf)]. 



Long-time accurate 2nd-order methods for Stokcs-Darcy systems 



25 







e u 




h 


BDF2 


AMB2 


BDF2 


AMB2 


BDF2 


AMB2 


1/16 


1.86e-005 


2.70e-005 


1.89e-005 


1.36c-004 


7.25e-003 


1.27c-001 


1/32 


1.71c-006 


2.17e-006 


1.67e-006 


1.02e-005 


1.88e-003 


2.59e-002 


1/64 


1.55e-007 


1.79e-007 


1.46c-007 


7.39c-007 


4.78e-004 


5.82e-003 


1/128 


1.38c-008 


1.51c-008 


1.28c-008 


5.47e-008 


1.21e-004 


1.38e-003 


favg 


3.46 


3.60 


3.51 


3.76 


1.97 


2.17 



Table 6.4 

Relative error and order of accuracy with respect to the spatial grid size h for Example 2 at 
t = 1 and with At = h. 



In this long time numerical experiment, we set the terminal time T = 100, and 
h = 1/64. We choose At = ^ for BDF2 and At = ^ for AMB2. The 
relative errors are plotted in Figures (|6.1[) - ([6.3[) . It is clear that although the errors 
grow initially, they remain bounded for all time. Moreover, the second-order in time 
accuracy for the velocity and the hydraulic head are also evident even in this onerous 
long-time experiment. The long-time accuracy for the pressure seems to be first-order 
in time, in agreement with our uniform in time error estimates. However, this is in 
contrast to the short-time second-order in time accuracy for p as recorded in Table 
1631 



Error of Hydraulic Head for BDF 



Error of Hydraulic h 



Fig. 6.1. Relative error for <f> in Example 3 for BDF2 (left) and AMB2 (right) up to t = 100 
for h = 1/64. 



Error of Velocity for BDF 



Error of Velocity for AMB 




FlG. 6.2. Same information as for Figure fcTTI but for Uf 
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Error of Pressure for BDF 



Error of Pressure for AMB 




Fig. 6.3. Same information as for Figure [6Jl but for p. 







e u 




h 


BDF2 


AMB2 


BDF2 


AMB2 


BDF2 


AMB2 


1/16 


2.05c-003 


2.95e-002 


1.49e-003 


1.72e-003 


4.88e-002 


1.70e-001 


1/32 


4.36e-004 


7.76e-003 


4.18e-004 


4.26e-004 


1.40e-002 


4.32e-002 


1/64 


9.84e-005 


1.99c-003 


1.09c-004 


1.07e-004 


3.64e-003 


l.lOc-002 


1/128 


2.32c-005 


5.05e-004 


2.75e-005 


2.68e-005 


9.29e-004 


2.79e-003 


Tavg 


2.15 


1.96 


1.92 


2.00 


1.91 


1.98 



Table 6.5 

Relative error and order of accuracy with respect to the spatial grid size h for Example 3 at 
t = 1 and with At = h. 



7. Concluding remarks. We proposed and investigated two long-time accu- 
rate and efficient numerical methods for coupled Stokcs-Darcy systems. The first is 
a combination of the second-order backward differentiation formula and the second- 
order Gear extrapolation method. The second is a combination of the second-order 
Adams-Moulton and Adams-Bashforth methods. Our algorithms are special cases of 
the implicit-explicit (IMEX) schemes. The interfacial term that requires communi- 
cation between the porous media and conduit, i.e., between the Stokes and Darcy 
components of the model, is treated explicitly in our algorithms so that only two de- 
coupled problems (one Stokes and one Darcy) are solved at each time step. Therefore 
these schemes can be implemented very efficiently and, in particular, legacy codes can 
be used for each component. 

We have shown that our schemes are unconditionally stable and long-time stable 
in the sense that solutions remain uniformly bounded in time. The uniform bound 
in time of the solution leads to uniform in time error estimates. This is a highly de- 
sirable feature because the physically interesting phenomena of contaminant seques- 
tration and release usually occur over a very long time scale and one would like to 
have faithful numerical results over such time scales. Spatial discretization is effected 
using standard finite element methods. Time-uniform error estimates for the Darcy 
hydraulic head and the Stokes velocity and pressure for the fully-discrete schemes are 
also presented. These estimates are illustrated by numerical examples. The methods 
proposed can be also utilized to approximate steady-state solutions in case the prob- 
lem data are time independent. All these features suggest that the two methods have 
strong potential in real applications. 
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On the other hand, there is still room for improvement. One could design even 
higher-order numerical methods. A third-order method was proposed in [10) with- 
out analysis. We are currently developing third-order unconditionally stable schemes 
based on the Adams-Moulton-Bashforth approach. It is also desirable to use different 
and adaptive time-steps for the two regions involved due to the disparate time-scales 
in the two regions that one sees in practical situations; see, e.g., (34J [35] . Also, mortar 
element method can be naturally adopted and may be useful to efficiently handle the 
different spatial scales in the two subdomains; see, e.g., [30] . 

So far, all methods deal with confined (saturated) karst aquifers. Most aquifers 
are unconfined and hence different methodologies involving either two-phase flows 
or free boundaries must be considered. Models for unconfined karst aquifers are 
inherently nonlinear. Mathematical investigation of unconfined karst aquifers is still 
in its infancy and deserves much needed attention. 

Last but not the least, the application of these methods to the quantification of 
uncertainty in flow and contaminant transport would be of great interest in real appli- 
cations that feature uncertainty in both the conduit geometry and matrix hydraulic 
conductivity. 
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